#Libraries
library(raster)
library(sp)
library(rgdal)
library(terra)
library(rgeos)
library(sf)
library(RStoolbox)
library(ggplot2)
library(gbm)
library(dismo)
library(rJava)
library(rmaxent)
library(rasterVis)
library(viridis)
library(tidyverse)
library(grid)
library(gridExtra)
library(magick)
library(plotly)
library(ecospat)
library(FactoMineR)
library(factoextra)

Population geography

Populations (n = 174) were chosen to span the entire distribution of Common tansy throughout the state of Minnesota. Common tansy is much more common in the northern half of the state, which is reflected in our sampling design. In the southern portion of the state, we attempted to collect plant tissue from as many reported populations as possible.

#Read in coordinates
pops <- read.csv("data/tansy_pop_info_mwsamples.csv")
xy <- pops[,5:6] #lat/lon of each of the populations in the dataset
coordinates(xy) <- ~ lon + lat
proj4string(xy) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")

# plot(xy, xlim = c(-93, -92), ylim = c(43, 50), pch = 19, cex = 0.1)
# textplot(xy@coords[,1], xy@coords[,2], words = rownames(tansy_tab), cex = 0.25, new = FALSE)

pops$geo_grp <- as.factor(pops$geo_grp)
pops$ECS_fac2 <- as.factor(pops$ECS_fac2)
pops_df <- pops

coordinates(pops) <- ~ lon + lat
proj4string(pops) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")

pops_df$finestructure_grp <- as.factor(pops_df$finestructure_grp)
pops_df$fs_clust_4 <- as.factor(pops_df$fs_clust_4)

#cluster 1 pops
clust1 <- which(pops_df$cs.K2.1>=0.5)

xy_1 <- xy[clust1]

#cluster 2 pops
clust2 <- which(pops_df$cs.K2.1<=0.5)

xy_2 <- xy[clust2]

Sample Eco-Regions

The MN Department of Natural Resources has divided the state using the Ecological Classification System using the National Hierarchical Framework of Ecological Units (ECS website). We assigned each population to a geographical grouping as indicated by the colors in the map. Groupings cluster populations that were similar in geography and ecology.

midwest <- map_data("county", "minnesota")
midwest <- midwest[midwest$long > -100 & midwest$long < -89.5 &
                         midwest$lat > 41 & midwest$lat < 52,]

MNsites <- ggplot() + geom_path(data = midwest, mapping = aes(x = long, y = lat, group = group), color = "gray") +
  geom_point(data = pops_df, mapping = aes(x = lon, y = lat, col = ECS_SUBSEC, group = name)) +
  scale_color_manual(values = ECS_cols) +
  coord_fixed(ratio = 1.5) +
  theme_bw()

ggplotly(MNsites)

Tansy spatially-explicit Genomic Structure

We used a conStruct analysis of populations to determine if there was spatially-explicit structure to genomic variation (GBS dataset of 3071 polymorphic loci). The analysis for K=2 again highlighted two geographically partitioned groupings. Group 1 (red) found mostly in the northeastern portion of the state, particularly St. Louis Co., along the North Shore, and into North Central MN. Group 2 has two distinct clusters - one in the Northwest portion of the state and the other in a portion of the forests of Lake Co. in the Northeast. Southern populations tend to be an equal sample of the two major groupings.

struct2 <- magick::image_read("conStruct/spatial/spK2_250K_structure.plot.chain_1.pdf")
#struct2_rot <- image_rotate(struct2, 270)
pie2 <- magick::image_read("conStruct/spatial/spK2_250K_pie.map.chain_1.pdf")

conStructK2 <- c(struct2, pie2)
image_append(image_scale(conStructK2, "x500"))

Environmental Variables

We downloaded and processed a broad range of environmental variables that may affect tansy population survival, growth, and dispersal for use in niche analyses. More detailed information about where data was obtained and how it was processed can be found in our GitHub

Niche Overlap analyses

Soil and Climate Niche

Hill-Smith analysis on mixed data types- aggregate environment

Ordination plot

#create background points for total climate
bkgrd <- spsample(MN_poly, 1000, type='stratified', iter = 1000)
projection(bkgrd) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")

names <- vector(length = length(xy))
names[clust1] <- "clust1"
names[clust2] <- "clust2"
IDs <- c(names, rep("bkgrd", length(bkgrd)))

MN_pts <- as_tibble(rbind(xy@coords, bkgrd@coords))
MN_pts$ID <- IDs

land_use_class <- raster::subset(MN_env_stack, 23)
env_varbs_analyze <- stack(MN_plant_growth_env_3_stack, land_use_class)

site_envs <- raster::extract(env_varbs_analyze, MN_pts[,1:2])
site_envs_1 <- as_tibble(site_envs)
site_envs_1$IDs <- IDs
site_envs_2 <- na.omit(site_envs_1)
site_envs_2$drainage <- as.factor(site_envs_2$drainage)
site_envs_2$land_class <- as.factor(site_envs_2$land_class)

pca_clim <- dudi.pca(site_envs_2[,1:3], scannf = F, nf=2) 
ecospat.plot.contrib(contrib = pca_clim$co, eigen = pca_clim$eig)

mix_env <- dudi.hillsmith(site_envs_2[,1:9], scannf = F, nf = 5)


ecospat.plot.contrib(contrib = mix_env$co[,1:2], eigen = mix_env$eig)

ecospat.plot.contrib(contrib = mix_env$co[,2:3], eigen = mix_env$eig)

ecospat.plot.contrib(contrib = mix_env$co[,c(1,3)], eigen = mix_env$eig)

Mixed CA scores

mix.scores.global <- mix_env$li

mix.scores.cl1 <- suprow(mix_env, site_envs_2[which(site_envs_2$IDs=="clust1"), 1:9])$li

mix.scores.cl2 <- suprow(mix_env, site_envs_2[which(site_envs_2$IDs=="clust2"), 1:9])$li

mix.scores.mn <- suprow(mix_env, site_envs_2[which(site_envs_2$IDs=="bkgrd"), 1:9])$li


mix.scores.tot <- as_tibble(cbind(site_envs_2$IDs, mix.scores.global[,1:2]))
colnames(mix.scores.tot) <- c("IDs", "Axis1", "Axis2")

biplot_env <- ggplot(data = mix.scores.tot, aes(x = Axis1, y = Axis2, color = IDs, shape = IDs)) +
  geom_point() +
  geom_point(data = subset(mix.scores.tot, IDs == "clust1"), aes(x = Axis1, y = Axis2, color = IDs, shape = IDs)) +
  geom_point(data = subset(mix.scores.tot, IDs == "clust2"), aes(x = Axis1, y = Axis2, color = IDs, shape = IDs)) +
  scale_shape_manual(values = c("clust1" = 19, "clust2" = 19, "bkgrd" = 21)) +
  scale_color_manual(values = c("clust1" = "green", "clust2" = "red", "bkgrd" = "gray80")) +
  theme_bw()

biplot_env

Testing for niche differentiation

Calculate and plot niche space

Axis 1 v. 2
mix.grid.clim.clust1 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,1:2],
                                          glob1=mix.scores.mn[,1:2],
                                          sp=mix.scores.cl1[,1:2], R=100,
                                          th.sp = 0,
                                          th.env = 0,
                                          extend.extent = c(-1,1, -1, 1))

mix.grid.clim.clust2 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,1:2],
                                          glob1=mix.scores.mn[,1:2],
                                          sp=mix.scores.cl2[,1:2], R=100,
                                          th.sp = 0,
                                          th.env = 0,
                                          extend.extent = c(-1,1, -1, 1))

ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2, 
                       title = "Niche overlap", name.axis1 = "Axis1", name.axis2 = "Axis2",
                                          col.unf = "#882255", col.exp = "#88CCEE", 
                                          col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
                                          transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,1:2], mix.scores.cl2[,1:2], mix.scores.mn[,1:2], mix.scores.mn[,1:2])

{cairo_pdf(filename = "figures/enm_niche_overlap_axes_1_2.pdf", height = 6, width = 8)
  ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2, 
                       title = "Niche overlap", name.axis1 = "PC1", name.axis2 = "PC2",
                                          col.unf = "#882255", col.exp = "#88CCEE", 
                                          col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
                                          transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,1:2], mix.scores.cl2[,1:2], mix.scores.mn[,1:2], mix.scores.mn[,1:2])
  }
dev.off()
## quartz_off_screen 
##                 2
Axis 1 v. 3
mix.grid.clim.clust1 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,c(1,3)],
                                          glob1=mix.scores.mn[,c(1,3)],
                                          sp=mix.scores.cl1[,c(1,3)], R=100,
                                          th.sp = 0,
                                          th.env = 0,
                                          extend.extent = c(-1,1, -1, 1))

mix.grid.clim.clust2 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,c(1,3)],
                                          glob1=mix.scores.mn[,c(1,3)],
                                          sp=mix.scores.cl2[,c(1,3)], R=100,
                                          th.sp = 0,
                                          th.env = 0,
                                          extend.extent = c(-1,1, -1, 1))

ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2, 
                       title = "Niche overlap", name.axis1 = "PC1", name.axis2 = "PC3",
                                          col.unf = "#882255", col.exp = "#88CCEE", 
                                          col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
                                          transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,c(1,3)], mix.scores.cl2[,c(1,3)], mix.scores.mn[,c(1,3)], mix.scores.mn[,c(1,3)])

{cairo_pdf(filename = "figures/enm_niche_overlap_axes_1_3.pdf", height = 6, width = 8)
 ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2, 
                       title = "Niche overlap", name.axis1 = "PC1", name.axis2 = "PC3",
                                          col.unf = "#882255", col.exp = "#88CCEE", 
                                          col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
                                          transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,c(1,3)], mix.scores.cl2[,c(1,3)], mix.scores.mn[,c(1,3)], mix.scores.mn[,c(1,3)])
  }
dev.off()
## quartz_off_screen 
##                 2
Axis 2 v. 3
mix.grid.clim.clust1 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,2:3],
                                          glob1=mix.scores.mn[,2:3],
                                          sp=mix.scores.cl1[,2:3], R=100,
                                          th.sp = 0,
                                          th.env = 0,
                                          extend.extent = c(-1,1, -1, 1))

mix.grid.clim.clust2 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,2:3],
                                          glob1=mix.scores.mn[,2:3],
                                          sp=mix.scores.cl2[,2:3], R=100,
                                          th.sp = 0,
                                          th.env = 0,
                                          extend.extent = c(-1,1, -1, 1))

ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2, 
                       title = "Niche overlap", name.axis1 = "PC2", name.axis2 = "PC3",
                                          col.unf = "#882255", col.exp = "#88CCEE", 
                                          col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
                                          transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,2:3], mix.scores.cl2[,2:3], mix.scores.mn[,2:3], mix.scores.mn[,2:3])

{cairo_pdf(filename = "figures/enm_niche_overlap_axes_2_3.pdf", height = 6, width = 8)
 ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2, 
                       title = "Niche overlap", name.axis1 = "PC2", name.axis2 = "PC3",
                                          col.unf = "#882255", col.exp = "#88CCEE", 
                                          col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
                                          transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,2:3], mix.scores.cl2[,2:3], mix.scores.mn[,2:3], mix.scores.mn[,2:3])
  }
dev.off()
## quartz_off_screen 
##                 2

Tests for overlap equivalency and similarity

#Overlap

D.overlap <- ecospat.niche.overlap(mix.grid.clim.clust1, mix.grid.clim.clust2, cor = TRUE)$D
paste(sprintf("The niche overlap is %f", D.overlap))
## [1] "The niche overlap is 0.421473"
#Niche equivalency

eq.test <- ecospat.niche.equivalency.test(mix.grid.clim.clust1, mix.grid.clim.clust2, rep = 1000,
                                          intersection = 0.1,
                                          overlap.alternative = "lower",
                                          expansion.alternative = "higher",
                                          stability.alternative = "lower",
                                          unfilling.alternative = "lower")

ecospat.plot.overlap.test(eq.test, "D", "Equivalency")

ecospat.plot.overlap.test(eq.test, "expansion", "Equivalency")

ecospat.plot.overlap.test(eq.test, "stability", "Equivalency")

ecospat.plot.overlap.test(eq.test, "unfilling", "Equivalency")

sim.test <- ecospat.niche.similarity.test(mix.grid.clim.clust1, mix.grid.clim.clust2, rep=1000,
                                          intersection = 0.1,
                                          overlap.alternative = "lower",
                                          expansion.alternative = "higher",
                                          stability.alternative = "lower",
                                          unfilling.alternative = "lower",
                                          rand.type = 2)

ecospat.plot.overlap.test(sim.test, "D", "Similarity")

ecospat.plot.overlap.test(sim.test, "expansion", "Similarity")

ecospat.plot.overlap.test(sim.test, "stability", "Similarity")

ecospat.plot.overlap.test(sim.test, "unfilling", "Similarity")

Climatic Niche

PCA - climate

PCA plot

#create background points for total climate
bkgrd <- spsample(MN_poly, 1000, type='stratified', iter = 1000)
projection(bkgrd) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")

IDs <- c(rep("site", length(xy)), rep("bkgrd", length(bkgrd)))

MN_pts <- as_tibble(rbind(xy@coords, bkgrd@coords))
MN_pts$ID <- IDs

site_envs <- raster::extract(MN_plant_growth_clim_3_stack, MN_pts[,1:2])

pca_env <- dudi.pca(site_envs, scannf = F, nf=2)
ecospat.plot.contrib(contrib = pca_env$co, eigen = pca_env$eig)

PCA scores

scores.global <- pca_env$li

scores.cl1 <- suprow(pca_env, site_envs[clust1, 1:3])$li

scores.cl2 <- suprow(pca_env, site_envs[clust2, 1:3])$li

scores.mn <- suprow(pca_env, site_envs[which(MN_pts$ID=="bkgrd"), 1:3])$li

names <- vector(length = length(xy))
names[clust1] <- "clust1"
names[clust2] <- "clust2"

IDs_spec <- c(names, rep("bkgrd", length(bkgrd)))
scores.tot <- as_tibble(cbind(IDs_spec, scores.global))

biplot_env <- ggplot(data = scores.tot, aes(x = Axis1, y = Axis2, color = IDs_spec, shape = IDs_spec)) +
  geom_point() +
  geom_point(data = subset(scores.tot, IDs_spec == "clust1"), aes(x = Axis1, y = Axis2, color = IDs_spec, shape = IDs_spec)) +
  geom_point(data = subset(scores.tot, IDs_spec == "clust2"), aes(x = Axis1, y = Axis2, color = IDs_spec, shape = IDs_spec)) +
  scale_shape_manual(values = c("clust1" = 19, "clust2" = 19, "bkgrd" = 21)) +
  scale_color_manual(values = c("clust1" = "green", "clust2" = "red", "bkgrd" = "gray80")) +
  theme_bw()

biplot_env

Testing for niche differentiation

Calculate and plot niche space

grid.clim.clust1 <- ecospat.grid.clim.dyn(glob=scores.global,
                                          glob1=scores.mn,
                                          sp=scores.cl1, R=100,
                                          th.sp = 0)

grid.clim.clust2 <- ecospat.grid.clim.dyn(glob=scores.global,
                                          glob1=scores.mn,
                                          sp=scores.cl2, R=100,
                                          th.sp = 0)

ecospat.plot.niche.dyn(grid.clim.clust1, grid.clim.clust2, quant = 0.25, interest = 2, title = "Niche overlap", name.axis1 = "PC1", name.axis2 = "PC2")

ecospat.shift.centroids(scores.cl1, scores.cl2, scores.mn, scores.mn)

Tests for overlap equivalency and similarity

#Overlap

D.overlap <- ecospat.niche.overlap(grid.clim.clust1, grid.clim.clust2, cor = TRUE)$D
paste(sprintf("The niche overlap is %f", D.overlap))
## [1] "The niche overlap is 0.443794"
#Niche equivalency

eq.test <- ecospat.niche.equivalency.test(grid.clim.clust1, grid.clim.clust2, rep = 1000,
                                          intersection = 0.1,
                                          overlap.alternative = "lower",
                                          expansion.alternative = "higher",
                                          stability.alternative = "lower",
                                          unfilling.alternative = "lower")

ecospat.plot.overlap.test(eq.test, "D", "Equivalency")

sim.test <- ecospat.niche.similarity.test(grid.clim.clust1, grid.clim.clust2, rep=1000,
                                          intersection = 0.1,
                                          overlap.alternative = "lower",
                                          expansion.alternative = "higher",
                                          stability.alternative = "lower",
                                          unfilling.alternative = "lower",
                                          rand.type = 2)

ecospat.plot.overlap.test(sim.test, "D", "Similarity")

ecospat.plot.overlap.test(sim.test, "expansion", "Similarity")

ecospat.plot.overlap.test(sim.test, "stability", "Similarity")

ecospat.plot.overlap.test(sim.test, "unfilling", "Similarity")